"""
Author: Pedro L. Magalhães
Email: pmlpm@posteo.de
Date: 2023
"""

#******************************************************************************
#******************************************************************************

# import libraries

import math

import pysoleng as pse

#******************************************************************************
#******************************************************************************
   
filename = 'tmy_54.027_9.803_2005_2020.csv'

(npdata, 
 names, 
 latitude, 
 longitude, 
 elevation) = pse.utils.read_jrc_tmy_csv(
     filename=filename)

# try to use the data

site = pse.timeplace.Location(
    coordinates=(latitude, longitude), 
    crs_string='epsg:4326')

site_ground_reflectance = 0.6

list_I_tilted_total = [] 
list_I_tilted_beam = [] 
list_I_tilted_sky_diffuse = [] 
list_I_tilted_albedo = []

for i, local_time in enumerate(npdata['timeUTC']):
    
    site_horizontal_total = npdata['Gh'][i]
    site_horizontal_diffuse = npdata['Gdh'][i]
    
    # find solar altitude and solar azimuth
    
    day_year = pse.solar.day_of_the_year(local_time, False)
    
    declination_angle = pse.solar.declination_angle_spencer1971(day_year)
    
    ast = pse.solar.lst_to_ast(local_time, site)
    
    solar_ha = pse.solar.hour_angle(ast)
    
    cos_zen_angle = pse.solar.cosine_zenith_angle(
        latitude*math.pi/180, 
        declination_angle*math.pi/180, 
        solar_ha)
    
    solar_az = pse.solar.solar_azimuth_angle(latitude*math.pi/180, 
                                             declination_angle*math.pi/180, 
                                             math.acos(cos_zen_angle), 
                                             solar_ha)
    
    solar_alt = pse.solar.solar_altitude_angle(latitude*math.pi/180, 
                                               declination_angle*math.pi/180, 
                                               solar_ha)
    
    # define the new surface

    surface = pse.timeplace.Plane(
        azimuth_zero_south_positive_west=solar_az, 
        slope_zero_horizontal_positive_south=math.pi/2-solar_alt)
    
    # compute the radiation on the tilted plane
    
    (I_tilted_total, 
     I_tilted_beam, 
     I_tilted_sky_diffuse, 
     I_tilted_albedo) = pse.models.tilted_plane_irradiance_isotropic(
         site, 
         surface,
         local_time, 
         site_ground_reflectance, 
         site_horizontal_total, 
         site_horizontal_diffuse)
         
    list_I_tilted_total.append(I_tilted_total)
    
    list_I_tilted_beam.append(I_tilted_beam)
 
    list_I_tilted_sky_diffuse.append(I_tilted_sky_diffuse)
 
    list_I_tilted_albedo.append(I_tilted_albedo) 
         
# compare the results

dni_pysoleng = sum(list_I_tilted_beam)*3600*2.7777778E-07
dni_jrcpvgis = sum(npdata['Gbn'])*3600*2.7777778E-07

print('Results for DNI:')
print('pysoleng (isotropic): '+str(dni_pysoleng)+' kWh/m2/y')
print('jrc-pvgis (?): '+str(dni_jrcpvgis)+' kWh/m2/y')
print('The error between them amounts to ' +
      str(abs(dni_pysoleng-dni_jrcpvgis))+' kWh/m2/y.')

#******************************************************************************
#******************************************************************************

# using dataframe

from pandas import DataFrame
import numpy as np
import datetime as dt
import pysoleng as pse
import math

class SolarDataFrame(DataFrame):
    
    KW_DAY_YEAR = 'doty'
    KW_DECLINATION = 'declination'
    KW_AST = 'ast'
    KW_SOLAR_HOUR_ANGLE = 'sol_ha'
    KW_COS_ZEN_ANGLE = 'cos_zen'
    KW_SOLAR_AZIMUTH = 'sol_az'
    KW_SOLAR_ALTITUDE = 'sol_alt'
    
    def __init__(self, filename: str):
        "A class for solar data for a given point on Earth."
        
        # TODO: allow users to provide the data manually
        
        # TODO: allow users to specify which correlations/models to use
        
        (npdata, 
         names, 
         latitude, 
         longitude, 
         elevation) = pse.utils.read_jrc_tmy_csv(
             filename=filename)
        
        #**********************************************************************
        
        # location(s)
        
        # TODO: variable location via GeoDataFrame
        
        site = pse.timeplace.Location(
            coordinates=(latitude, longitude), 
            crs_string='epsg:4326',
            elevation=elevation)
        
        # add columns to numpy array
        
        #**********************************************************************
        
        # solar position columns
        
        sol_pos_array = np.ndarray(
            shape=(len(npdata),), 
            dtype=[
                *npdata.dtype.descr, # add tmy data
                (self.KW_DAY_YEAR,float),
                (self.KW_DECLINATION,float),
                (self.KW_AST,dt.datetime),
                (self.KW_SOLAR_HOUR_ANGLE,float),
                (self.KW_COS_ZEN_ANGLE,float),
                (self.KW_SOLAR_AZIMUTH,float),
                (self.KW_SOLAR_ALTITUDE,float)
                ]
            )
        
        # TMY data
        
        for dtype_tuple in npdata.dtype.descr:
            
            sol_pos_array[dtype_tuple[0]] = npdata[dtype_tuple[0]]
        
        # day of the year
        
        sol_pos_array[self.KW_DAY_YEAR] = tuple( 
            pse.solar.day_of_the_year(local_time, False) 
            for local_time in npdata['timeUTC']
            )
        
        # declination
        
        sol_pos_array[self.KW_DECLINATION] = tuple( 
            pse.solar.declination_angle_spencer1971(day_year)
            for day_year in sol_pos_array[self.KW_DAY_YEAR]
            )
            
        # apparent solar time
        
        sol_pos_array[self.KW_AST] = tuple( 
            pse.solar.lst_to_ast(local_time, site)
            for local_time in npdata['timeUTC']
            )
        
        # solar hour angle
        
        sol_pos_array[self.KW_SOLAR_HOUR_ANGLE] = tuple( 
            pse.solar.hour_angle(ast)
            for ast in sol_pos_array[self.KW_AST]
            )
        
        # cosine of the zenith angle
        
        sol_pos_array[self.KW_COS_ZEN_ANGLE] = tuple( 
            pse.solar.cosine_zenith_angle(
                latitude*math.pi/180, 
                declination_angle*math.pi/180, 
                sol_ha)
            for sol_ha, declination_angle in zip(
                    sol_pos_array[self.KW_SOLAR_HOUR_ANGLE],
                    sol_pos_array[self.KW_DECLINATION]
                    )
            )
        
        # solar azimuth
        
        sol_pos_array[self.KW_SOLAR_AZIMUTH] = tuple( 
            pse.solar.solar_azimuth_angle(latitude*math.pi/180, 
                                          declination_angle*math.pi/180, 
                                          math.acos(cos_zen_angle), 
                                          sol_ha)
            for sol_ha, declination_angle, cos_zen_angle in zip(
                    sol_pos_array[self.KW_SOLAR_HOUR_ANGLE],
                    sol_pos_array[self.KW_DECLINATION],
                    sol_pos_array[self.KW_COS_ZEN_ANGLE]
                    )
            )
        
        # solar altitude
        
        sol_pos_array[self.KW_SOLAR_ALTITUDE] = tuple( 
            pse.solar.solar_altitude_angle(latitude*math.pi/180, 
                                           declination_angle*math.pi/180, 
                                           sol_ha)
            for sol_ha, declination_angle in zip(
                    sol_pos_array[self.KW_SOLAR_HOUR_ANGLE],
                    sol_pos_array[self.KW_DECLINATION]
                    )
            )
        
        DataFrame.__init__(self, 
                           data=sol_pos_array, 
                           columns=sol_pos_array.dtype.names)
        
        self.site = site
        
        # self.sol_pos_array
        
        # TODO: associate units with each column
        
        # x = np.array([(1,2),(3,4),(5,6)],dtype=[('a','<i4'),('b',float)])
        
        # DataFrame.__init__(self, data=x, columns=x.dtype.names)
        
        #**********************************************************************
        
        

class CollectorSurfaceDataFrame(SolarDataFrame):
    
    KW_SURFACE = 'surface'
    KW_INCIDENCE_ANGLE = 'incidence'
    KW_IRRAD_PLANE_TOTAL = 'g_plane_total'
    KW_IRRAD_PLANE_BEAM = 'g_plane_beam'
    KW_IRRAD_PLANE_SKY_DIFFUSE = 'g_plane_sky_diffuse'
    KW_IRRAD_PLANE_ALBEDO = 'g_plane_albedo'
    
    def __init__(self, filename: str):
        
        # TODO: multiple surfaces
        
        site_ground_reflectance = 0.6
        
        # surface(s)
        
        surf_array = np.ndarray(
            shape=(len(npdata),), 
            dtype=[
                ('surface',pse.timeplace.Plane),
                ('incidence',float),
                ('g_plane_total',float),
                ('g_plane_beam',float),
                ('g_plane_sky_diffuse',float),
                ('g_plane_albedo',float)
                ]
            )
                
        # surface
        
        # TODO: fixed surface 
        
        surf_array['surface'] = tuple( 
            pse.timeplace.Plane(
                azimuth_zero_south_positive_west=solar_az, 
                slope_zero_horizontal_positive_south=math.pi/2-solar_alt
                )
            for solar_az, solar_alt in zip(
                    sol_pos_array[self.KW_SOLAR_AZIMUTH],
                    sol_pos_array[self.KW_SOLAR_ALTITUDE]
                    )
            )
        
        # TODO: incidence angle (irrelevant with tracking)
        
        # sol_pos_array['incidence'] = tuple( 
        #     pse.solar.cosine_incidence_angle(declination_angle, 
        #                                      latitude, 
        #                                      plane_slope_horizontal, 
        #                                      plane_azimuth_south, 
        #                                      hour_angle)
        #     for solar_ha, declination_angle in zip(
        #             sol_pos_array['sol_ha'],
        #             sol_pos_array['declination'],
        #             sol_pos_array['sol_alt'],
        #             sol_pos_array['sol_alt']
        #             )
        #     )
        
        # irradiation/irradiance on the surface
        
        for t, time in enumerate(npdata['timeUTC']):
            
            (surf_array['g_plane_total'][t], 
             surf_array['g_plane_beam'][t], 
             surf_array['g_plane_sky_diffuse'][t], 
             surf_array['g_plane_albedo'][t]
             ) = pse.models.tilted_plane_irradiance_isotropic(
                  self.site, 
                  surf_array['surface'][t],
                  time, 
                  site_ground_reflectance, 
                  npdata['Gh'][t], 
                  npdata['Gdh'][t]
                  )
                 
        # self.npdata = npdata
        # self.sol_pos_array = sol_pos_array
        # self.surf_array = surf_array
        
        # TODO: concatenate ndarrays into one to create the dataframe
                 
        # npdata = np.concatenate((npdata, sol_pos_array, surf_array), axis=1)
        
        #self.npdata = npdata
        
        # create dataframe
        
        # DataFrame.__init__(self, data=npdata, columns=names)
            
        
pse.classes        
        
filename = 'tmy_54.027_9.803_2005_2020.csv'
sdf2 = SolarDataFrame(filename=filename)

#******************************************************************************
#******************************************************************************